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Abstract 

Spectral nultigrid nethodo are demonstrated to be a ccmpecitit^e technique 
for solving the transonic potential flow equation. The spectral 
discretization, the relsKation scheme, and the nultigrid techniques are 
described In detail. Significant departures from current approaches are first 
illustrated on several linear problems. The principal applications and 
examples, however, ore for compressible potential flow. Tnese exatiples 
Include the relatively challenging case of supercritical flow over a lifting 
airfoil. 
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I* 

Interest in tra^tiionic acrcaynamics endurec because rsoet ccnrserclal and 
military aircraft operate prodcmlnantly in the transonic reglno. Tlic design 
nnd analysis of transonic uingo and related configurations have bean carried 
out largely within the fra^aework of the transonic small perturbation equation 
and the full potential equation. Apart from their relative simplicity the 
popularity of these flow models is due to their adequate representation of 
flow features of practical importance. For Instance, the pressure rise across 
an Iscntroplc shock in these models is sufficiently accurate for normal Jiach 
numbers ahead of the shock less than 1.3. Naturally, if other design 

considerations produce strong shocks and/or complex vortical flow, then 
recourse to the Euler equations is appropriate. Indeed, Euler solutions to 
transonic flow problems have attracted serious attention as of late, and they 
surely gain inc.eai.±ng popularity as they become more competitive with 
potential solutions. However, for many configurations of engineering 
lnter.-st, potential flow predictions with asymptotically first-order weak 
vlscous-invlacld interaction give solutions of more than adequate accuracy 
tlj. i;hen strong shocks and/or vortlclty are of dominant importance in the 
flow field, weak vlscous-lnvlseld interaction is no longer an adequate 
model. Implementation of strong interaction models is relatively crude at 
this time, and until substantial improvements have been made, the potential 
formulation will retain the most favorable accuracy-to-cost ratio for a wide 
range of practical transonic flow problems. 

The main difficulty in the numerical solution of the steady transonic 
flow problem has been the mixed elliptic-hyperbolic nature allowing for the 
presence of discontinuities. The initial breakthrough in overcoming this 
difficulty was made only in the early 1970's by Ifuman and Cole [2] who 
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Port„r.atlon conation. fol.onln, t„d„ .rca,=t„ro„r.n thorn have heen can. 

^cvclopccnta 1 „ the conpntnt.on of trnnaon.c flo-oo. ^o carve. Icctaron of 

Bnllhana and d.caon ,4, orocent a detailed revlcn of the.c dc.cloprentc 
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factorization schemes, AFl and av '> t ^ 4 

, Ulch, applied to the transonic EJdall 
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the approximate factorization scheme AF3 developed by Baker [8] 
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dlscnsse, the b.sle ptloelple el the appterlnate factotlrntlon nehe.es fet the 
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Doria and South [11 J. Another fast method Is the inultlgrld technique, first 
npplled b. south and Ht.ndt (.2, to the ttansonle snail petturbatlon equation 



with SLOR as a baul.c iteratlva scheua. Recently, Jaussoa (13] developed the 
ciultlcrld procedure to accelerate convercence of the full potential solution 
by an ADI aethod. Despite the CKistence of quite a few efficient nathodo of 
potential colution, controlled coznparisonc are lacking. 

The conputer tine required to obtain nunerical solutions for tvo- 
dinenolonal potential flows ic now eo saall that there is practically no 
Incentive for developing nore efficient schenes. However, for threc- 
dloensional flows existing methods are still so costly that a substantially 
more efficient solution algorithm would have great practical importance. 
Unlike the two-dimensional case, computer storage is a crucial consideration 
in weighing the efficiency of a scheme. Paeudospectral methods have 
demonstrated their capacity for producing equivalent accuracy with far fewer 
grid points than standard second-order or even fourth-order methods, not only 
for smooth flows but also, more recently, for the Euler equations [14]. The 
first paeudospectral two-dimensional potential flow solutions were obtained by 
Streer.t (15], who established that equivalent solutions ware in fact obtained 
for potential flows with far fewer grid points than required by standard 
methods. However, his- solution technique was clearly in need of acceleration, 
particularly for supercritical flows. In this paper we describe an 
acceleration technique, based on the spectral nultlgrld methods developed by 
Zang, ct al. (16], (17], that has significantly improved the rate of 
convergence of the pseudospectral discretization of the full potential 
equation. In fact, the spectral multlgrld scheme is so efficient that the 

preliminary version described here Is highly ccmFcCltlve with the finite 
difference schemes* 

Since the application of spectral methods to compressible flov/s is still 
e fairly novel approach, most readers are likely to bo unfamiliar with either 
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the practical details of opactral tsathoda or the nuances of r.as;erlcol cethods 
for coapreoslble flows. Horeover, spectre! nultip.rld nethodo thenoelvee are 
otill in the formative stage. The promialng nature of the preoent results 
warranto a reasonably complete and self-contained deocriptlon of the nunierical 
Taethod. 

We begin by describing a means of implementing pocudoopcctral differenti- 
ation, which, although asymptotically inefficient, is nonetheless prefcroble 
for problems on moderately-sized grids. This is followed by descriptions of 
the essential features of spectral multlgrid methods and of the relaxation 
schemes. These methods are then illustrated on several linear problems. An 
explanation of the potential flow problem and its pseudospectral approximation 
is given next. Finally we report on the performance of the spectral multlgrid 
method on both aubcritlcal and supercritical potential flows. 


II. Spectral Kathode Uaiag Katrls Kultlplco 
Ihe Fast Fourier Transform (FFT) has usually been cited as a key element 
in the efficiency and hence the implementation of spectral methods. In the 
pseudospectral sort of calculations discrete Fourier methods are commonly used 
in the evaluation of derivatives. However, under some circumstances it is 
actually faster to use conventional matrix-vector multiplications for this 
purpose than to resort to transform techniques. An obvious requirement is 
that the problem be of moderate size. There arc many significant engineering 
applications which meet this requirement. The transonic tlow application, 
which Is the main thrust of this paper, is one such example. Even in circum- 
stances which most favor transform techniques — on grids with 2^' points 
the matrix-multiply approach (using nothing but Fortran) has proven 
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to be sicniflccntly faster than tha transfonn matiiod (enploylns ascembly 
language FFT'o) . Prccisa coraparloona will be given below. 

In a pceudoopcctral Eethcd the fundamental reprsoentatlon of the solution 
is In physical space. The quantities t+.lch are stored are the values of the 
function u(x) at special collocation points Xj. Derivatives, however, are 

evaluated spectrally. Tne values of the function are passed through a 
suitable discrete transform to produce the repreaentatlon of the function in 
transform (wavenumber) space. The actual differentiation takes place in 
wavenumber space. Then an inverse transform is applied to yield the 
pseudoopectral approximation to tha derivatives of the function at the 
collocation points. Let U denote the vector of values of the function at 

the collocation points. Then the approximation to the derivative at these 
points may be written 

0 u, (1) 

where 

0 - C”^ D C, (2) 


with C representing the discrete transform and D representing 
i^tlon In wavenumber space. 

The most well-known pseudospectral method is baaed upon Fourier series. 
Let tha Interval of Interest be [0,2n] and use the collocation points 

2irl 

“ M J (3) 


Then 
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(4) 





6 


PACE ^ 

OF POOR QUAL.Tlf 



k 

I 


M 

2 

M 

2 




1 



( 5 ) 


(c-M 


jk 


2irllf,i 


c 


M 


(6) 


The Fourier series differentiation matrix may be conatructed by the matrix 
multiplies implied by Eq. (2). Alternatively, one may simply use the explicit 
formula given in Eqs. (8) and (9) of [16] for the elements of 0. 

Once the matrix 0 has been constructed, the cost of evaluating a 
derivative by the matrix vector product OU is O(M^). The transform 
technique reduces this to 0(M £n M) . However, two transforms are required 
and the constant in the 0(H £n H) factor is larger than the one In the 
O(M^) case. 

Chebychev pseudospectral methods have been the most widely used ones for 
non-periodic boundary conditions. The standard Interval is [-1,1] and the 
collocation points arc 


tM 

X. “ C09 
J N 




(7) 


Then 




" Hi 7 N 

k-J 


k,j “ 0,1,«**,K , 
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■■ r ’ — 
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where 


Moreover, 


where 


and 
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j « 0 or N 
otherwise 



I > k+1 and 1 = k+1 (nod 2) 
otherwise 



j » 0 
otherwise 



cos 



(9) 


( 10 ) 


(ID 

( 12 ) 


An explicit formula is available in Eqs. (49) and (50) of [16] for this 
Chcbyshev differentiation matrix* 


III* Spectral Kultigrld PundciaGntalo 
Overview of Hultl^rid Aly^orlthms 

The problems of Interest here are scalar partial differential boundary 
value problems. The PT)S can be written in the general form 

L (u) ” f , (13) 

where u(x,y) is the unknown function, f(x,y) is some source term, 
and L is a partial differential operator which might be nonlinear in the 
unknown u. The correaponding diocrote problem will be written 
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L(;i) « F (14) 

In obvious notation. 

Hultlgrld solution schones for Eq. (l/>) involve conbinlng relaxation 
sweeps for that equation with relaxation sweeps for related problems on 
coarser grids. Let V denote an approximation to U. The essential property 
for the relaxation scheme is that it preferentially damp Che high-frequency 
components of the error V - U« Then after a small number of relaxations the 
error will have so little high-frequency content that it can be approjslnated 
well on a coarser grid. Solutions on the coarser grid are relatively 
Inexpensive to obtain, especially if this strategy is applied recursively by 
using still coarser grids as needed. 

Let us consider just the interplay between two grids. The fine-grid 
problem is written 

L^(U^) " F^. (15) 

The shift tc the coarse grid occurs after the fine-grid approximation has 

been sufficiently smoothed by the relaxation process, l.e., after the high- 
frequency content of the error heo been sufficiently reduced. The 

related coarse-grld problem is 


L^(U^) «= F^, 


(16) 


where 

r[f^ - L^(V^)] + L*^(RV^), (17) 

The restriction operator R interpolates a function from the fine grid to the 
coarse grid. The coarse-grld operator and solution are denoted by L^ and 
respectively. After an adequate approximation to the conrsc-grid 
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problen has been obtained, the flnc-srid approjtimatlcn 1 g corr 


ected vifi 


+ P(V^ - liV^), 


TJia prolongation onerator P interpolates a function fro= the coarse grid to 


the fine grid. 


The choice of the coarsc-grid problem Ir based upon rewriting Eq. (1 


5) as 


<■ [f^ - L^(V^)] + L^(V^). 


The tern in brackets is the fine-grid residual. Since it has been presuned to 
be smooth, its coarse-grld approximation is clearly 

Equations (16) ond (17) th.n (oHou by teplaolns the rot,.lnl„s flnc-stld 
quantities with appropriate coarse-grid ones. 


The quantity 


- RV^ 


is the coarse-grid correction. Equations (16) to (IB) are equivalent to 


L^(RV^ + W'^) - L*=(RV^ « 


+ PZ*^, 


where is the 


3ppro:'.lmation to W^. For linear problems Eq. (22) rcdi 
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■ F^. 


(24) 


This overview has, of course, bser. based upon the paper by Krannt [18J , 
albeit In notation popularized by Hackbusch [19J. The particular choices of 
the Interpolation and coarse-grld operators used in the present spectral 
multlgrld work are described In the following sub-sections. This description 
Is given for one-diKcnslonal probleas. The extension to higher dlncnslons Is 
obvious. These details are followed by a discussion of the relaxation 
schemes • 


Interpolation Operators 

The spectral multigrid interpolation operators vhich were proposed in 
[16] for periodic coordinates amount to trigonometric interpolation: given a 

function on a coarse grid (v»lth points), compute the discrete Fourier 
coefficients and then use the resulting discrete Fourier series to construct 
the interpolated function on the fine grid (with points)* This nay be 

accomplished by performing two FFT's. An explicit representation of the 
prolongation operator is 

(25) 

(26) 


If 


- 1 


P B . 


f C 


«i « — ■ ■■ 1 

2 ^ ‘ 


which sums to yield 


p m — Sf^ - —1 
Jk M M 

etc 


where 



r. ; 


u 


S(r) 


:ln(vfrM^)cot<irr) 


Or PCCil v^-i. 


coo(TrrM^) 


J:\JT 


r integer 


otherwise 


( 27 ) 


The corresponding restriction operator is esnentlally the adjoint c. 


thlo: 


n .. i_ sfJ_ _ 1 l 

Jk % I., 

3 C f 


■). 


rs) 


Interpolation for non-periodic coordinates employs Chebyshev scries in an 
analogous fashion. Tlie prolongation operator is 


jk 


i - _i 

c. 


c, N 

k c 


vlV 

cos ^ cos 

f c 


where Cj^ Is detlned by Eq. (9) with N - N . This s 


sums to 


( 2 ?; 




where 


(30) 


Q(r) 


r integer 


V 4 - V 4 cos(irrN^) + l/j cos(|^(N^.(-l ) ) slnf-^] csc(— ) 


■Y”-) otherwise 


(Jl) 


We will have occasion to use two distinct restriction operators. One is 
sometimes used in forming the coarse-grid operator and Is obtained bv applying 
Chebyshev restriction In the obvious fashion. It will ho denoted bv r('i) 

Qnd It In given by 


R 


(o) ^2 r- 

jk — 


Ir)K 


* IT 

c '-f 


(32) 
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where Cj^ Is defined by Eq. (9) with N » 11- and 


N 

Va+ — 


V4 + V2 co3(^(N^+l))sin[ — 


r Integer 


othen/ine 


The other Is used for Interpolation, is denoted by and Is defined 

the adjoint requlreoent: 




c ‘ f "c "f 


where c. is defined by Eq. (9) with I? « N . 

•» r* 


Coarse-Grld Operator 


A typical term in the class of problems considered here is 


^ [a(u,x) 


The discrete operator which represents its fine-grid pscudospectral 
approximation is 

- 0 A 0, /ofiv 


diagonal matrix 
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Many nultigrid Invaatlgctors* c*g., [1?], [20], and [21], have advocated 
choosing the coarsc-grid operator go that 

« RL^P. (3o) 

Both the Fourier and the Chebyshav flrst«-derivative operators, defined by Eqs* 
(2) - (12) , satisfy 

0^ « RO^P, (39) 

where R * is chosen in the Chebyshev case* However, Fq* (38) itself is 

not satisfied if the coarse-grld analog of Eq. (36) is used to define 
except In tha trivial case for which a(u,x) is a constant* On the other 
hand, much of the efficiency of the pseudospectral method is lost if Eq. (38) 
is used to define the coarse-grid operator. Some compromises were suggested 
in [17]* The most satisfactory one seems to be using Eq* (36) but with the 
restricted values of a(uj,Xj) in place of tha pointwlae values* Tlie 
Chebyshev restrictions should be performed with R^^^* 

Boundary Conditions 

In the applications? that follow, three types of boundary conditions 
appear: periodic, Dlrlchlet, and Neumann. Periodic boundary conditions are 

automatically satisfied by the use of Fourier series, ^"ully-perlodic problems 
contain some subtleties that are discussed in [17]. 

Dlrlchlet boundary conditions are handled effortlessly. The vector of 
unknowns should include the values at the boundary points in their natural 
locations. (This has th(> side, effect of facilitating the programming of the 
Chebyshev Interpolation.) On the fine grid the desired boundary values are 
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inserted Into the appro[>riate loc£tioa»3 and these values nre not codified 
during the relaxation* On the coarser grids the epproprinte boundary vcIugs 
are the ones tihich fall out of the restriction process. 

Neunann boundary conditions are a bit touchier. Wc have enforced then by 
incorporating the Ncunann boundary condition into the discrete operator. 
Suppose that there is a Ueumann boundary condition at x » -1. In the 
evaluation of a tern euch as appears in Eq. (35) , the first stage is the 
computation of du/dx at all the collocation points. In general this value 
will not match the desired boundary value. Hie boundary condition Is enforced 
by resetting the value of du/dx at x ” —1 to the desired value before 
proceeding with the multiplication by a(u«x) and then the final 
differentiation. This produces the desired boundary condition in the 
converged solution. This approach has the advantage of ensuring that the 
boundary condition appears in the discrete operator with a consistent scaling. 
A much less effective alternative is to replace the differential equation at 
X “ -1 with the condition that du/dx is the prescribed boundary value. The 
disadvantage of this approach is that this boundary equation is far out of 
scale with the rest of the operator. This alternative has in fact beea tried 
on some of our test problems and it has resulted in a substantial 
deterioration of the convergence rate. 


IV. Felaxatloa Schescs 

The crucial property that a relaxation scheme should possess for use in a 
nultigrld algorithm is that it damp effectively the high-frequency components 
of the error. It need not be especially effective in the low-frequency range, 
GO long ns it docs not: iimplify any coTponents . For spectral multigrld methods 
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™ additional ro,olro,.ont orl.e. £t<» th. „l„bol ntturo of tt,o opprotit^tlon: 
tho foot ov.loot.oo Of dorivotlvo. dcoodo thot tho roteotlon bo oiooltoocoo. 
rotbor tboo ootoooolvo, o.g., jooobj.,, it.plo=o„tod offtoiootly. 

whereas Gauss-Seidel's cannot. 

o. d-teratlve ochemes that macto these requireraents is based upon 

approximate factorisation techniques C5]. n.ese nethoda are especially 

attractive because they have been employed in come of the most successful 

finite difference solutions to the delicate transonic potential flov problem 

t7], [13] . Moreover, the latter work demonstrated their effectiveness in tho 

multigrid context, albeit for a purely finite difference approximation. A 

review of the computational transonlcs literature suggests that the most 

fruitful interpretation of approximate factor! ». 

dce ractorlzcitlon schemes for this nixed 

elliptlc-hypevbolf. probl™, i. .orr..po„dl„p tte.-d.po„d.„c 

partial dltf.r.ntlal .,„atl„„. Thia is th, .pp.o.ch .n.t will bs ,.k.„ below. 

An alternative end perhaps eoro tr.dltlenal interpretation tor linear, 
elliptic problems la in terms of preconditioning. roe relaxation schcite 
p oposed in II 7J for a spectral multigrld method for aach problems was 
iaterpreted as aa laeo.plete decompositloa serving a. a preeoaditloaiag for 
Aiehardsoa-a iteration. A briaf deseriptioa of thia scheme is iaeladed her. 

aiaee it will serve as a co.p.rlsoa for the approximate f.ctorlratloa method 
on one of the linear test problems. 


Ric hardson Iteration with I ncomplete LU necompo.sl Mnn 

A preconditioned Richardson iteration for solving Kq. 
expressed as 

V ■«- V + 10 h"^ [F - L(V)] , 


(14) can be 


(40) 
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whore n is the preconditioning r-ntrix and tj io the roloKotlon parameter. 
Tlie matrix H should be chosen r.o that it io an approximate inverse to L, 
but is easily invertible. Ihc vernlon recommended In [17] for linear problems 
is obtained by first constructing the matrix Hpj, vhlch represents a standard 
second-order finite difference approximation to L (sec Eq. (13)) and then 
performing an incomplete LU decomposition of I?pp. Details are provided in 
[17J along with a prescription for choosing the relaxation parameter u bo 
that the high-frequency error components are damped preferentially. 


Approximate Factorization 

For this discussion it is convenient to rewrite Eq. (14) as 


where of course. 


M(U) - 0, 


M(U) = L(U) - F. 


( 41 ) 


(42) 


Next, view U not as the solution to Eq. (41), but rather as the steady-state 
solution to the evolution equation 

a7“M(u). 

This is surely sensible If L(u) Is elliptic for then Eq. (43) represents 
the spatial discretisation of a parabolic problem. Scml-lmpllclt time- 
stepping procedures are desirable for such problems because of the severe 
explicit time-step limitations. (This la especially acute for pseudospectral 
discretizations employing Chebyshev series because of the very small spacing 
between the collocation points near the boundary.) Tl,o simplest practical 
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ttea diocretisatlor. of Eq, (43) is POOS QUALITY 



„(n+l) _ y(n) 

At - + J(U^"^)(U 

where 




(n+1) . 

(44) 




(45) 


and a superscript refers to a tine level. Let 




a - 1- 



V 

and 


(46) 


r. u^n+l) _ y(n)^ 


(47) 

- 

and then rewrite Eq. (44) as 



j \ 

[al - J(U^"^)]au^"^ » 

). 

(4n) 


where I denotes the Identity matrix. 




This motivates the relaxation scheme 



f* 

V + V + uAV, 


(49) 

t 

where AV Is the solution to 



?• # 
i 

f- 

f 

f 

i. 

[al - J(V))AV - M(V). 


(50) 

1 

t 

V 

In many cases the Jacobian J(V) can be split Into 

the sum of two 

operators 

f 

dj{(V) and dy(V), each Involving derivatives In 

only the one 

coordinate 
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direction Indicated by th« subocrlpt. Approsimte factorisation cotbodo 
encompass various approxltiatlono to the left-hand cide of Eq. (48). The moot 
straightforward of those la 

** foil - Jy(V)JAV « aM(V), ( 51 ) 

in combination with Eq. (49). This is Just the Douglas-Gunn version of ADI 
[22]. It is commonly referred to as AFl for the transonic problem [5]. For 
second-order spatial dlecretizatlons the terra [al - J (V)] leads to a set of 
tridiagonal systems, one for each value of y. The second left-hand aide 

factor produces another set of trldlagonal systens. For pseudospectral 

discretizations, however, these systeraa are full; hence, Eq. (51) la ctlll 

relatively expensive to Invert. A corapromlse analogous to the one Invoked In 

the Incomplete LU decomposition preconditioning Is to replace J and J 

X y 

with their second-order finite difference analogs > denoted by and H.,, 

respectively: 


[al - H (V)] [al -- H (V)]AV « aM(V) . 

X y 


(52) 


The approxinate factorization scheme consists of Eqs* (49) and (52). For 
purely finite difference approximations some analytical results are available 
for selecting optimal values for the parameters a and w' (9). Ko similar 
results arc yet available for the present application. By analogy with the 
finite difference case we have chosen (o to be of order unity and have 
selected a sequence of a's In a range by the rule 


o, 

h\x ^ 


k-i 
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shore K denotoe the nothot of dlotlnet o'n. The cfcofcoe of o^ end o,^ 
sore based In port on eetfeete, of the elsenvnlno mn;;o of the dloereL 
oporatero end fn <n„ch creoter) port bp trim end error, rortenatelp. the dn 
echene Ip not very oensitivc to these paratacters. 

for olnsle-ijrld ooletlons to the onbcrlticm pctentiel flo, problen the 

pseudoopectral AFl achene baaed on E". (L-i^ h-„ . . 

hciO proven satiofcctorj’’ [15]. 

intensive verb on finite difforenee nethod. for oepererltlcnl potontlnl floe 
hn. indlceted the necooelty to been their opprotlnote foetorltetlon echeneo on 


2 

3^C 

do3t 


M(U), 


(54) 


"here e 1. n physlenl vnrleble directed mens the etreenllne. One echen, 
nhtch nodole thin behavior lo referred to no AFl (SJ. A pecudoopectral a'2 
vnrlent Is described In (13,. since schene. of the AFl type nodel hyperbolic 
eqnetlono they ere relatively Ineffective at donpins hlsh-freqncncy error 
coeponents. Indeed. In the pseodospoctral elnglc-grld Inpleientotlono (15, 
for enpercrltleal flop, on Iterative .trategy Involving both API and AF, 
found to be uore effective than AFl alone. (By Itself, of c.nrse, AFl van 
divergent.) Ihl. mil be referred to belov as the AFl/AFl scheae. 


V. nsrsarlcal ncoulto for Linear Problcaa 
We chose a scries of test problems to bridge the gap hetveen the spectral 
.eltlgrld method, described In ,17, and tho.c re,nlred for the potential flow 
ptoblem. The first step was to change the relegation schc.e from 
preconditioned Richardson fcetatlon to appconlcato Coctotieotlon. Tho 
boundary conditions veto loft es Dlrlehlct In both coordlo.te directions. The 



nest phaoG involved chiftlns Co periodic boundary eonditiono in one 
direction. In the final atar»^ the gecaetry too altered frc 3 a rcctcnsle to on 
annuluo t>lth on inaei: radial boundary condition of Heunenn rather th?»?i 
Dirlchlet type. Tliis lest problca is about ac close ao one con ccr^c to the 
potential flow problen uithin the constraint of linearity. 

The nultlgrid codes used a casimja of 4 levels. These are labelled by 
the index k, where k«2, 3, or 5. The grid cn level k contalnn 
either 2^' (Fourier) or 2^ 1 (Chobyshev) collocation points in a coordinate 
direction (including boundary points). Two different schedules were used; 
they were referred to as schedules B and D in il7]. For schedule B the 
problem was first solved on level 2; then that soluticc was interpolated to 
level 3 as the initial guaos for a cultigrid iteration Involving levels 2 and 
3; then the converged level 3 solution was interpolated to level 4 as its 
initial guaca, and so on until level 5. For nclicditlc D the multigrid process 
simply began on level 5. In both cases the initial guf^ss consisted of random 
numbers chosen from (0,1)^ ensuring that all error components were present 
Initially. Both schedules were run in a fined code with 6 relaxations (2 
passes turough a 3 parameteir sequence) before restriction to a coarser grid. 
A coarse-grid solution was deoned acceptable for prolongation to a fine-grid 
whenever its PUS residual dropped below O.IZ of the last residual on the finer 
grid. All of these linear runa employed the correction scheme, i.e., Eq. (24) 
rather than Eq. (16) was solved on the coarser levels. The variable 
coefficients and the right-hand sides for the coarsc-grid problems vrere 
filtered in the manner described In [ 17 ] • 

The specific measure used was the equivalent ancothlng rate. In some 
preliminary cclculationa the average time required for a single fine-grid 
relaxation was determined. For an actual tiultigrid calculation let and 
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T 2 be tha EKS rccidualo after the first and last flr.c-srld relcrationo, 

respectively and let t be the total CPU ti=e. Then the equivalent eaoothins 

rate rao taken an 


1 



(55) 


Rectangular Chebyahev - Chebvohev Prohlen 

The problen class is the same one examined la [17]: 


— Ini J. 3 r 3ui 

8x1“ drJ + « f. 


(56) 


on (-1,1) X (-1,1) with Dlrlchlet boundary conditions with 


a(y.,y) «» 1 + e 


COS TO TT(x+y) 


(57) 


and f(::,y) and the boundary data chosen so that the solution is 


^(x,y) “ sln(m^xx + r/4) oln(n^iry + n/4), (58) 

The properties of three test cases are listed In Table I. The parameters used 
in the approximate factorization scheme are given in Table II. 

The performance of the preconditioned Richardson (PR) and the approximate 
factorization (AT) methods Is shown In Table III. Tbe PR method Is sbout 
twice as fast as AF on these problems. But recall that the PR scheme hr-s been 
highly. tuned (especially for problem 1), whereas the ,\F scheme was subjected 



22 




.. only n nnnll »n„nl o.- „lnl an. nrror tunln,. „n ,onbn tHn „ n^nenn 
™ul. innoan ,„nnly Innn n„o onpnnlnnnlnilnn, non no nn„il„„ .nnlynln. „n 

content with ectebllshing Ito wrkabllity In thlo nultlerid contest. 

imon dorlvntlvo conlnntlons nto perfomod via FFT'p, tho tlnn coqnlpnd on 

CDC Cyber 175 for e 5 pelejiation (Includlnp both tlie resldnel 

evalnerlon end ,ecrorlrn„» prp,„, ,, 

only nbonr 51 ol rbe rornl rlne In rb.n, celcnlnrlnnn ,pe„r 
polnrlnp berneen leveln. On ...„pn there „ere 4 to 5 rele.eelnn. ,er every 

interpolation. A ccsparlcon between the tr~nnf« j 

ecween the tranoforn and natris-multlply nethods 

ol dlllerentlntlon In provided In y.ble XV. Only level 5 ,n 33 e 33 prld, 

one. one pnln by veins FFr.. P„rtbe„,ere. sine, .nnt m tbe verb tnbee pleee 
on level. 2 ro 4. the t.t.l mnnlns ti.e 1 , 1 ... (by 10 . 2 „„ 

P > verslona. Bear in mind that assembly language FFT'o were performed 

on grids ideal for the FFT (powers o<^ Tho 

KP ers o- .). TTie matrix ra.ltlplles were coded in 

Fortrnn. In tbe petent.nl llev nppli.nrlen It 1. edventnpeevn te verb en vere 
oooernl ,rlde. .^v. the „tr.r-™lt,piy eltern.t.ve 1. highly tenpetltlve. 
1.0 ndvnntege evght to ertend te even Inrger grid, en vector prece.eers. 


Table I. G,aractcrlotlca of the Rectangalar 
Cliebychov-Chebyshcv Teot Probleaa 


Problem No. 

e 

'"u 


— 

1 

0.00 

1 

1 


2 

0.20 

2 

2 


3 

1.00 

5 

5 
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TTablc lie Parastttoro o£ tha A? Schss::® for the 

Chsli3?CrtGV‘’*dlct>yoliCV Pfolsl f 


Level 


“h 

(0 

2 

1 

6 

1.4 

3 

8 

75 

1.2 

4 

80 

1000 

1.1 

5 

600 

8000 

1.0 


Table III- Eqelvalcnt Sasothins T^tcn on the 

Eectcngalar Chcbyahev-Chcbysliev Problesas 

Problem No. 

FR 

AF 

1 

.26 

.43 

2 

.58 

♦ 78 

3 

.78 

.92 


Table IV. Hcoidual Braluation Tice for the AF Schese on the 
Rectaneulor Chcbyaliev-Chebyabev Teot Problesuj 


Level 

Transfonn Method 
Differentiation 

Matrln-Kaltlply 

Differentiation 

3 

.0204 

.0083 

4 

.0622 

.0390 

5 

.214 

.248 
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Rectangular C hebyshcv-Fourter Prohlpn 

Thl« ptoble. IS also dsscnbsj by a,. (j„_ ^ 

»lth Dirlchlst boundary conditions In a and periodicity In y. ibc 
coefficient 


a<x,y) «= 1 + e e 


cos nig{TTx+y) 


(59) 


and the rest of the problem fits the solution 


u(x,y) - 3ln(n^iix + r/4) slnCm^r cos y + r/4). 


(60) 


The propcrtlc, of three teat ca.e. are Hated In Table V and the AF parameter, 
are supplied In Table VI. 

Table VH gives the results. There la evidently nothing to be gained 
here by .orbing np to the finest level by first solving the coarser level 
problens. The present eonblnatlon of the eoarse-grid operator and the AF 
paraneters .onld not pernit a solution to be obtained for a highly oscillatory 
proble. such as the previnns sob-seetlon's problen 3. Pete that the 

equivalent snoothlng rates on the present problons 2 and 3 are eo.parable to 
those for the previous problem 2. 


Table V. Charactcrlotlcn of the Kectangular 
aiebyshcv-Fourler Teat Problces 


Problem No. 

c 


•"a 

1 

0.00 

1 

1 

2 

0.10 

1 

1 

3 

0.20 

2 

2 
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Tabic VI. ParascCcro o£ tbc I'S Jtebcsc for tbc 

Ecctcacular Chebyahcv-Fourler rccbleca 


Level 

“l 

% 

(U 

2 

0.5 

6 

l.O 

3 

2,0 

75 

1.0 

A 

10. 0 

1000 

1.0 

5 

100.0 

8000 

1.0 


Table VII* Equival^^nt Ss:ootLinp, ilatco on 
Uectcngular Cticbyohev—Fourlcr 

tbe 

Problc!2f3 

Problem No. AF/B 

AF/D 

1 .77 

.75 

2 -78 

.79 

3 .82 

.76 


•I 


Annular Chchyshev-Fourlof Problem 

The differential equation for this last linear example is 





[7 


« -^1 


f 


on (1,5) X (0,2x) with 


cos(rn (r+0)) 
a 


(61) 


V 


a(r,0) = 1 + c c 


( 62 ) 
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The radial boundary conditions are T^lrichlct at r =* 5 but TTsunann at r » !• 
Periodicity In czinuth is enforced. Tables VIII and IX present the test case 
and AF paraneters, respectively. 

The results are available in Table X. Tnese are the least impressive 
smoothing rates of the linear test problems* Ncuvnann boundary conditions are 
usually more troublesome than Dirlchlet ones* The global approx Imatlcn 
urderlylng the spectral methods makes them especially difficult to enforce* 


Table VIII* Character is tics of the Annular 
ChcbyGhcv-Fouricr Test Pr ob lens 


Problem No. 

a 


“a 

1 

0.00 

1 

1 


0*10 

1 

1 

3 

0.20 

2 

2 


Table DC. 

Parameters of the 
Annular Giebyshcv- 

AF Scheme for the 
-Fourier Problems 

Level 


a (1) 

h 

2 

5 

40 2*0 

3 

10 

600 1.4 

4 

100 

6000 1.0 

5 

iOOO 

10000 1.0 
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Tcbla S. Giuivalcnt £?3cothlng Eateo on tho 
/anolar Gsebyshev-Fouricr Probloss 


Problem No. 

AF/B 

AF/D 

1 

.82 

.87 

2 

.81 

.87 

3 

.87 

.86 


Potential FI<kj Past an Airfoil 

Tl.a probl™ considered Is ihct of co.presolblc potenrlsl flow post o t»o- 
dlmeoslosol airfoil. Be „odel this rith the fell potential e, nation, applylns 
boendary conditions at the actual airfoil surface. In this port a numerically 
generated conformal mapping I23J m „ transfom, the nlrfoll onto the 

unit circle. The form of the transformation betue.n the complex physical 
plane (the t-plane) and the complex computational plane (the o-pl.ne) Is 

,(l-n) 




r% 


(64) 


coefficients A., and are generated numerically so that the 

knopu relation, betoeen the surface tangent angles nnd arc lengths of the 

airfoil shape are satisfied. The trailing edge of the airfoil Is located nt 

0-1 In the computational plane. The Schuart-Chrlstoffel factor In the 

transformation allows the smooth mapping of a finite-angle trailing edge. Tor 
further details on this particular mapping see Jameson I23J. Tlia Inner 
portion of a 16 x 48 grid is shoun in Figure 1. 
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In tbs coaputational plane, with c « Rc , the potential equation 
beccass 

Is (""ll!* lot! !)■'>. («> 

where <J Is the velocity potential end p is the density, given by the 
isentroplc relation 


1 

„ fi Y~I w2, 2 . 2 ,xlY-l 

p = [1 - — M„(q^ + qg - D] 


( 66 ) 


the ratio of specific heats is denoted by Y, the Ibch number at infinity is 
denoted by M^, and the velocity components in the physical (r,0) plane are 


with 



i n 

n 

H 3R 


I 

CD 

a 

RH 1 

H “ 

ill 

doj * 


The boundary conditions at the surface and in the farficld are 


li 

3R 


at R - 1 


•<* R cos 0 + E tan ^ [ v^l - tan o] 


as R 


(67) 

(63) 


(69) 


(70) 

(71) 


TTie first term in the farfleld boundary condition describes the uniform 
freestream flow. The remaining term Is the first-order lifting term; It is 
derived in [2A1. The quantity E is known as the circulation. It is 



ft 


r 

f-r 
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deten^lned by the Kutte condition, vhlch states that the physical velocity at 

the tralllne edge nuct be finite. Since H - 0 at the trallinE edce. the 
Rutta condition reduces to 
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3-> ^ 


at <T = 1« 


(72) 




The singularity of the potential In the farfleld poaca dlfficultieo 
(especially for spectral methods) that are best handled by computing in terms 
of the reduced potential G, vhich is defined by 


G " « - (R + |) cos 0 - E tan"^[/l - tan ©] 
and is assumed to be periodic In 0. It follows that G satisfies 


(73) 


along with 


aR-* 8 o'r 30^ “ 


(74) 


9G 

G 0 


at R « 1 
as R 


(75) 

(76) 


and the Kutta condition. 

The spectral method employs a Fourier series representation in 0. 
Constant grid spacing in G corresponds to a convenient dense spacing In the 
physical plane at the leading and trailing edges. The domain in R (with a 
large, but finite outer cutoff) is mapped onto the standard Chcbyshev domain 
[-1.1] by an analytical stretching transformation that clusters the 
collocation points near the airfoil surface. The stretching is so severe that 
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the retJo ,£ che l.rseet-t.-..,Ue.t redial laioi-vel. le „,er .000 lor tie 

8tid oho .0 moor portion 1. lllnotretod 1. rip„ro .. Too tr.nolornatlon le 

Incorporated Into the operetor oMoh reprceents differentiation In tile R 
direction. 

"ecplte It. nonllnenrlty the potential flon proUet, ronotn. fa.ny 

stralchtforward co lon» as thp finr. -r 

the floT-. is everyvihere nubsonlc. The real 

difficulty of tho problen an... uh.n the flo» fonts a oupor.onlc bobble on 

the .l.toll, Tbe potential oou.tlon 1 . then of nixed elllptlo-hyporbollc typo 

ond .dnlt. OO.R dl.contlnn.tlo,. Both oonproo.lon end 

enpan.lon .boob. .Ill .ppe., ^ 

bias is Introduced into the equation -fn 

quatlon in the supersonic region. The nost 

expedient technique for dealing wltb -hd. d„ 

witn IG to use the artificial dennlty 

approach of Hafez, et al« F2S1 j , 

[5]. The original artificial density is 


with 


P - P <5 p 


V - maxfo.l - 

JI 


(77) 


(78) 


I, th. local .b,oh nnnbor and Jp firot-otdor 
(ondlvldod, dlffotcnoo. 1 „ tbo pro.oot notb . b,„bor-otdot ottlflolol dcnolty 
fomuU tolatod to o tom dovolopod by J.oo.on 1.3) h.. boon mployod. 

Iho first p.ondo.poctr.l solution, to tho ccprosslblo potential fl.„ 
Ptobla, „oro Obtained by Sttaott I. 5], ,26, nslnp o slnglo-prld vor.lon of the 
cpproxlnat. feotorloetlon Itorotlv. sob... do.orlb.d In the fouttb sootlon. 
Por .ubotltloal flo.s tUls n.tbod ». elrondy b.pbly conpetltlyo mtb stetc- 
Of-tbe-ntt finite dlffotonoe notbods. P,r sup.rcrltlo.l flo.s. bo.ovor. tbo 
slnBlo-Ptld p.oudospeotr.1 eobooo ,„fro m.ff lol.nt. oyon «tb tbo uso of 



the AF2 extension of the approximate factorisation ochetae* Thio problems 
then, poses a useful applicatioiv of the spectral raultigrid approach and, as 
the results Indicate, a dranatic demonstration of its effectiveness. 


VII* Haoulto for roteatial Fltnf Faat ca /drfoll 

Tlic numerical cxamplea of this section have been chosen primarily to 
illustrate the effectiveness of the multigrld approximate factorization 
(MG/AF) solution scheme in comparison with the earlier single-grid approximate 
factorization (SG/AF) method [15] for solving the spectral equations for 
potential flow. A secondary Issue is the comparative quality of this spectral 
discretization and of wldely-uscd finite difference approxinatlons# A by- 
product of these examples Is soma practical guidelines for the multigrid 
algorithms. 

Three test problems suffice for a comprehensive treatment of the spectral 
multigrid efficiency and spectral discretization accuracy Issues: a 
Gubcrlcical lifting airfoil, a supercritical nonlifting airfoil, and a 
supercritical lifting airfoil. These have been listed In order of Increasing 
difficulty. Detailed comparison of the spectral SG/AF and liG/AF schemes will 
be provided for the first two examples. Extensive comparisons are also made 
for all three problems between the spectral MG/AF scheme and tv/o popular 
finite difference codes: TAIR [7], a slngle-grid/AF2 method and FL036 [13], a 
mul tlgrld/AF method. 

Some of the relevant issues have already been discussed in [15]. The 
most sensitive matttr is surely the weighing of the efficiency of two schemes 
(spectral and finite difference) with different accuracy and convergence 
properties. The reader is directed to [15] for a more detailed discussion 


than is provided here. 
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Three different Qvido have been used (with the coarser levels In 
parentheses): 16 x 32 (12 x 16 and 8 x S). 16 >= 48 (14 x 32, 12 x 16 and 
8 X 0) and 18 x 64 (16 x 48, 14 x 32. 12 x 16 and 8 x 8). Note that in 

passing to a coarser level the grid is typically reduced by less than a factor 
of 2 in each coordinate direction. -nus choice leads to a sienlflcant 
improvement over the standard griddlng for the spectral potential flow 

problem, especially in the supercritical regime where the solution has large 
high-frequency content. 

This problem has the added complication of a highly-stretched grid in the 
radial direction. This Is accounted for by changing the spectral 
differentiation matrices from C~^DC (see Eq. (2)) to 

0 = BC-1dc, 

where B is a diagonal raatrij: which contains the Jacobian of the 

transformation. A substantial improvement in the spectral rcultigrld algorithm 
results from defining the coarse-grid differentiation matrices directly by Ev. 
(39) rather than by the coarse-grid version of Eq. (79). In the absence of 
stretching these two definitions are equivalent. Equation (39) is easily and 
efficiently implemented with matrlx-raultlply techniques. 

Virtually all the spectral multigrid results Included here were obtained 
with the same fixed schedule: start on the finest grid, work down to the 

coarsest grid and then back up to the finest grid; on the way down there is 1 
sweep chough the (three) parameter sequence and on the way up there arc 2 
sweeps • 
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lubcrltical Llftintr Airfoil 

Tho flov past an KACA 0012 airfoil at angle of attack and a freeotream 
Mach nunber of 0.5 ulll serve as the first toot case. Tne airfoil produces a 
fairly large lift coefficient at these conditions and the surface preonure 
distribution shows a sharp suction peak near the leading edge. Since the 

local L^ch nuaber in this peak is nearly 1 , cotapresslbllity effects arc 
oubotcntlal* 

In order to denonstrate that the spectral solution on a relatively coarse 
grid captures all the essential details of the flow ve first compare It with 
an extremely accurate finite difference result. Tn Figure 2 la shown the 
surface pressure coefficient from a spectral solution using 16 points In the 
radial (R) direction, and 32 points In the azimuthal (0) direction; the 
symbols denote the solution at the collocation points. For comparison, the 
result from the finite difference code FL036 is shown as a solid line. The 
grid used In the benchmark finite difference calculation is so fine (64 x 384 
points) tnat the truncation error Is well below plotting accuracy. The 
spectral calculation seems to lack detail near the leading edge suction 
peak. However, since the spectral solution Is actually a continuous 

representation of the solution. It may be expanded In terms of its basis 
functions onto a much finer mesh. Such an expansion, shown in Figure 3 . 
reveals the hidden detail of the solution. The FL036 and expanded spectral 
results are identical to plotting accuracy. The spectral computation on this 
mesh yields a lift coefficient v^th truncation error less than 10 “^ Spectral 
solutions on a 16 x 32 are thus of more than adequate resolution and 

accuracy for subcrltlcal flows. 

The convergence histories for both the SG/AF and the MG/AF spectral 
schemes on this test case are displayed In Figures 4 and 5. The 


convergence 
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hlator-»ed hove been cupplled for both the icaKiuua residual (Fieure 4) and the 
error In circulation (Figure 5). They are plotted agalnot raachlne time on a 
CDC Cyber-175 cotaputcr. Although the rr^ltlgrid cede (henceforth referred to 
as •'MGAI'SP”) shows a substantial improvement over the single-grid appro::lm.ate 
factorisation code (”AFSP-) m ma..lmue residual convergence, the gain is even 
norc draraatlc from the lift convergence standpoint. Tnlo Is understandable 
since the lift is predominantly a low-frequency property of the solution. The 
single-grid spectral approximate factorisation schema was recognized to be 
weak In damping for long-wavelength error components [15]. 

The consensus In the computational transoulcs community appears to be 
that TAIR Is the fastest widoly-available finite difference code. A 
comparison of maximum residual versus machine time for TAIR and KGAFSF on the 
Bubcrltlcal test case Is shown in Figure 6. The two codes require nearly 
equivalent machine time with TAIR showing a better asymptotic convergence 
rate. However, the TAIR result was produced on a rather coarse (default) 
finite difference mesh of 30 x points. Compared with the surface pressure 
results from MGAFSP and FL036. the TAIR result Is significantly in error near 
the leading edge (Figure 7). This is Indeed truncation error, because TAIR 
results on a 60 x 297 mesh are more In agreement with those of HGAFSP and 
FL036. A further Indication of the somewhat large truncation error of the 

TAIR result Is that the predicted drag and lift coefficients are correct to 
only two decimal places (subcrltlcal potential flow yields identically zero 


drag) 


In Figure 8 are shown convergence histories from TAIR, FL036, and IfGAFSP 
on meshes which yield approximately equivalent accuracy; the surface pressure 
results are the same to plotting accuracy, the lift coefficient is converged 
In the third decimal place, and the predicted drag coefficient is less than 
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.001 • (Actually, the spectral result ia an order of magnitude core accurate 
than these Units, but the TAIR result barely siectc thera.) As can be seen 
from Figure A, the single-grid AFOP rccult uould fall in the vicinity of the 
FL036 and TAIR results in the present figure. 

Use of more Chan three grids in the spectral r.ultigrld code did not yield 
an improvement In effective convergence, since the interpolation overhead 
becaica a greater proportion of the total work. It would have been desirable 
to use a lowest grid coarser than 8 8 in the rnultigrid cycle. 
Unfortunately, due to the prcoence of the metric singularity at the trailing 
edge, coarser mssh results vere no oscillatory as to provide no useful long- 
wavelength Information. 

Supercritical Tlonllfting Airfoil 

The test is again the NACA 0012 but at « 0.8 and id.th zero angle of 
attack, l.e., a nonlifting condition. The surface pressure coefficient 
distribution as computed by the spectral method on an 18 x 64 grid is 
displayed in Figure 9. The shock at mid-chord Is relatively strong; the 
normal Mach number ahead of the shock is approximately 1.25. The shock is 
spread over several mesh spaces by the finite difference artificial viscosity 
used in the spectral calculation. Although this sliock is already far sharper 
than those produced by finite difference codes on a comparable grid, it ought 
to be possible to capture the shock In a still smaller region with a spectral 
method employing an artificial viscosity more suited to .the spectral 
discretization. 

The coi.vergence histories for the SG/AF scheme (combining AF2 and AFl) 
and the KG/AF scheme (using AFl alone) on a fine grid are shown in Figure 
10. The multigrid scheme obviously shows a much higher as>m:iptotic convergence 



rate. Kote that the slngle-srld achorae initially oocillatas with the maxinun 
rcoidual of order unity for a rather lengthy period. Thia lo indicative of 

the lack of high-frequency daaping in ths AF2 scheme. Ihc flcvfleld io being 

cstabliohcd in this period by the AF2 scheme; the plot of the history of the 

number of supersonic points in Figure 11 shows that the AF2 scheme establishes 

the shock position and the size of supersonic region nearly as fast as the 
nultlgrld scheme, albeit trtth some transient overshoot. This rapid 
establishment of the flcwfleld is at the expanse of high-frequency error, 
which is subsequently damped wlien the AF2/AF1 alternate cycling is begun. Tne 
nultlgrld algorithm, however, monotonically establishes the flowfleld and 
damps high-frequency errors in a far more efficient manner. 

Experience with all forms of transonic potential flow calculations has 
shown that convergence rates are quite sensitive to the order and amount of 
artificial viscosity: more artificial viscosity generally yields faster 

convergence, but at the expense of more Tridely smeared shocks. Multigrid 

schemes have been especially sensitive to these effects, and the present one 
la no exception. However, the large Improvement in efficiency offered by the 
multlgrld over the previous single-grid spectral scheme has allowed the use of 

much finer grids, offsetting the present, uncomfortably large artificial 
vlscoGlty . 

Supercritical Lifting Airfoil 

The lifting supercritical tent case was the NACA 0012 at » 0.75 and 

o 

a - 2 , which yields a section lift coefficient of nearly 0.6. A shock 
appears only on the upper surface for these conditions and Is rather strong 
for a potential calculation; the normal Mach number ahead of the shock Is 
about 1.36. Lifting supercritical test 


cases are especially difficult for 



opectrcl ccthois clnce tha colution will nlxfayjj have eignificant content In 
tha entire frequency spectrum; the chock populates the hlgheot frequencies of 
the grid end the lift ia predcaincntly on the ccale of the entire dosain. /n 
Iterative cchcrsa therefore cuot be able to dcap eri*or cosponanto across the 
spcctrun. Tha AF2/AF1 schene of [15] uas socevhat unreliobie for such 

problcns; so a coapuriuon vrlli not be oho^m between A?2/AF1 and the nultigrid 
schene. 

A history of the surface pressure coefficient le supplied in Figure 12. 
This dcKonotratCB the rapid convergence of the entire frequency cpectrua of 
tha solution. Pressure dlntributions arc cho’un after 0, 1, 4» and 9 cycles of 
the fired“cyclc algorithm; one cycle requires approximately 5 seconds of 
Cyber-175 titsa. The shock overshoot seen in tha 4-cyclc frame is a phenomenon 
associated with tha final positioning of the shock by the multigrid scheme. 
The finite dlfferance nultigrid cchc^Be exhibits clnilur behavior [13], 

All of ♦•he supercritical spectral nultigrid calculations shoxm thus far 
used a sequence of five rather than three grids, mostly due to the finer 
finest grid used for these casco. Scheduling within the fixed-cycle nultigrid 
algorithm was much the sane as for the subcritlcal cases: one or two passes 
through the tine-step sequence were made on each grid. Convergence for 
supercritical cases Is not always nonotonlc because adjustments in lift or 
shock position can Introduce high-frequency errors which nay require an extra 
cycle to damp. An adaptive cycle algorithm night be of benefit here provided 
that the "limit cycle" problem were avoided. 

Surface pressure distributions, both at the collocation points and 
spectrally expanded onto finer spacing, are shown in Figures 13 and 14 for 
grids of 16 X 48 and 18 x 64 points, respectively. As can be soon, the 
coarser-grld result predicts virtually the same shock position an the finer- 
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grid ccsputation; the lift coofficlcnto agree to IZ* These recults may be 
coGsparcd with those fraia the finite difference codec, TAIR and FL036, ohovm in 
Figurec 15 end 16, respectively. Tnn chock predicted by TAIR is for norc 
rounded end cleared than that of FL036, reflecting the coercer nsoh and larger 
artificial viscosity used in the fomer. The TAIR result shown is also only 
correct to one dec Inal place in lift as cenpared \?ith a fincr-grid result. 
Convergence histories for these four cases: spectral nultlgrld (16 x 48) 

and (18 x 64), TAIR (30 149), and FL036 (32 x 192) arc shown in Figure 17. 

The spectral results arc obviously handicapped in this comparloon by the 
necessity of such fine (for spectral nethods) neshes brought about by the use 
of the finite difference artificial viscosity form. Perhaps the purely 
spectral shock— capturing methods currently under development will pcrnlt sharp 

I 

shocks to be captured with still coarser meshes. 

VIII. Concluolons 

Spectral nultlgrld methods are still In their Infancy. Nevertheless, 
they have already exhibited the capacity to accelerate drastically iterative 
schemes for nonlinear, as well as linear, problems. Pough estimates of the 
asymptotic convergence rates indicate that the multlgrld procedure has led to 
an ireprovement over the single-grid spectral method of nearly a factor of 10 
for subcrltlcal cases; the improvement is considerably greater for 
supercritical situations. 

The worth of the spectral discretization Itself for compressible flows is 
now clear: equivalent solutions are Indeed obtained with far fewer gild 

points than are required for finite difference solutions. Since oubcrltical 
flows are smooth, the present results, showing both that the spectral method 
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convergence rate Is far better than cocond-erder and also that its absolute 

error level la lower than finite difference ones even on unreasonably coarse 

ends, are no ourFrlse. Undaniably, any shock discontinuity In supercritical 

flow should produce some dccradatlon In the formal accuracy of uhe spectral 

solution. Ronetheless, grid refinement studies demonstrate that the spectral 

solutions stabilize on far coarser grids than do finite difference 

solutions. Coupled with multlgrld solution techniques, spectral methods for 

steady compressible flows have reached the stage at which they are truly 

competitive with finite difference methods on problems of aerodynamic 
interest • 

Several aspects of this technique have to be Improved before spectral 
methods for compressible flows reach their full maturity. The preaent 
rela.ratlon schemes are just straightforward modifications of the ones used for 
finite difference methods. c^^ely relaxation schemes more tuned to the 
spectral discretization can and will be devised. There is also the clear need 
to develop more suitable forms of artificial viscosity for capturing shocks by 


spectral methods* 
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Figure 1. Grid about NACA 


0012 airfoil 
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Figure 2. Comparison of surface pressure 
coefficient distributions 
Spectral: 16 x 32 points 

FK036: 16 x 192 points 

NACA 0012, M - 0.5, a - 4° 



Figure 3. Surface pressure coefficient fron 
c:cpandcd spectral solution 
16 X 32 points 
IfACA 0012, II 0.5, a - 4° 
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Figure 4. Maximum re£3ldual vs. machine tine 
multigrld and single-grid schemes 
HACA 0012, K « 0.5, a = 4° 





Figure 5 


Krror in lift vs. machine time 
multigrid and single-grid ncheimes 
KACA 0012, H - 0.5, a - 4° 



lnf .,0 llrlL \' 


TAIR 
30 X 149 


MGAFSP 
16 X 32 
3 GRIDS 


8 . 12 . 16 . 

T (sec - CY175) 


- 2.0 


- 1.6 


- 1.2 


-f- TAIR 
PL036 




/ -2 .4 .6 

^ x/c 'V 


Figure 6. Maxlrcun residual vc. machine time 
NACA 0012, M » 0.5, a « 4° 


Figure 7. Cocparlson of surface prec£ 
coefficient dlstributlona 
TAIR; 30 x 149 points 

FL036: 32 x 192 points 

??ACA 0012, H « 0.5, a = 4° 
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Figure 8. Maxltaua residual vs. nachlne time 
NACA 0012, M « 0.5, a - 4° 



Figure 9. Surface proaourc coefficient 

opectral eolutien: 18 x 64 points 

I7ACA 0012, M 0.0, a ° 0 
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Figure 10. Haxlnum rcoldual vs. nachlne tine 
multigrid and single-grid echeoeo 
NAC* 0012, H - 0.8, a » 0 
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Figure 11. Humber of supersonic points 
VO* machine tima 
MCA 0012, M « 0*8, o M 0 
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a« Initial condition b. After 1 cycle 
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Fif^urc 12. Hintory of 8u 
sjpectral solu 
NACA 0012, M • 



d. After 9 cycles 


face prcosure coefficient diatrlbution 
Ion: 18 X 64 points 

0.75, a - 2” 
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c. Collocation points b. Eitpanded 

Figure 13. Surface preooure coefficient distribution 
spectral solution 16 x 48 points 
NACA 0012, K » 0.75, a « 2° 



a. Collocation points b. Expanded 


Figure 14- Surface pressure coefficient distribution spectral 
solution: 18 64 points 

KACA 0012, M » 0.75, a = 2° 




Figure 15« Surface pressure coefficient distribution Figure 16* 

TAXR: 30 x 149 points 

NACA 0012, M « 0.75, o - 2 ° 
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FL036: 32 x 192 points f. y: 

NACA 0012, M o 0.75, o - 2° n 1,; 

3 1-5 



53 


OF FCOn QuALiTV 



Figure 17 


Maximuia residual vs. machine time 
NACA 0012, H - 0.75, a 2° 








